Part of the Accelerator Readiness Review is looking at muons that leave the ring. To this end, we (Renee, Nathan, James, and myself) have a preliminary version of MDC-2 with full injection, unified fields, and “ghost” detectors in a Geant parallel world (so these sensitive volumes can be placed without overlapping physical structures). We can remove various magnets and fields and see how that affects muons escaping the beam.
# Load libraries
library(readGallery)
library(stringr)
library(dplyr)
library(ggplot2)
library(rgl)
library(readr)
library(purrr)
library(parallel)
The useDataProduct calls should all come together too.
useDataProduct('std::vector<gm2truth::TrackingActionArtRecord>')
useDataProduct('std::vector<gm2truth::RingTrackingPlaneArtRecord>')
useDataProduct('std::vector<gm2truth::GhostDetectorArtRecord>')
I have run preliminary MDC-2 (using Renee’s FCL files), making 10K events for various scenarios. Files are currently stored in /pnfs/GM2/scratch/users/lyon/arr_20170313 and directories within. Ghost detectors included are Renee’s cylinder just at the outer vacuum wall and encompossing the inflector as well as my ghost detector that is a cylindrical shell just on the inside of the world cube. Note that this code includes the magnet yoke steel.
Need a function to turn a /pnfs path into a xrootd url
xrootdify <- function(p) {
# We /pnfs/BLA --> root://fndca1.fnal.gov/pnfs/fnal.gov/usr/BLA
paste0('root://fndca1.fnal.gov/pnfs/fnal.gov/usr', str_replace(p, '/pnfs', ''))
}
# Let's look at only the 180* job series (so we just get 10,000 events per scenario)
system("ssh lyon@gm2gpvm04.fnal.gov 'ls /pnfs/GM2/scratch/users/lyon/arr_20170313/180*/*.root'", intern=T) %>% xrootdify() -> arrFiles
arrFiles
[1] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root"
[2] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089095_0/ARR_unified_noDipole_cyl.root"
[3] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089097_0/ARR_unified_noInflector_cyl.root"
[4] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089099_0/ARR_unified_noKickerQuads_cyl.root"
[5] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089123_0/ARR_unified_noKicker_cyl.root"
[6] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089132_0/ARR_unified_noQuads_cyl.root"
Let’s pluck out the scenario.
scenarios <- str_match(arrFiles, 'arr_20170313/.+_unified_(.+)_cyl')[,2]
scenarios
[1] "everything" "noDipole" "noInflector" "noKickerQuads" "noKicker" "noQuads"
What data are in these file? Well, product_sizes_dumper doesn’t seem to work across XRootD. That’s too bad.
Here’s a picture from Root.
Create a reader class (modify outside of this notebook) for TrackingActionArtRecord
readerClassSkel('gm2truth::TrackingActionArtRecord', writeFile = 'trackingActionReader.py')
trackingActionReaderFor documentation, here is the TrackingActionArtRecordReader reader class
readr::read_file('trackingActionReader.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class TrackingActionArtRecordReader(GalleryReaderBase):
def __init__(self, inputTag, evLimit):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID',
'parentTrackID', 'volumeUID', 'status', 'x', 'y', 'z', 'px', 'py', 'pz']
self.evLimit = evLimit
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))
def fill(self, ROOT, ev):
if ev.eventEntry() > self.evLimit:
return False
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2truth::TrackingActionArtRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2truth::TrackingActionArtRecord
for e in p: # Loop over elements and fill
# Let's only accept the muon
if e.trackID != 1 or e.parentTrackID != 0:
continue
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID, e.parentTrackID, e.volumeUID,
e.status, e.x, e.y, e.z, e.px, e.py, e.pz ])
return True
Make an instance
trackingActionReaderClass <- createReaderClass_from_file('trackingActionReader.py')$TrackingActionArtRecordReader
Load the data - let’s just look at one file since it takes a long time
taReader <- trackingActionReaderClass(artInputTag("artg4"), 20000)
getGalleryData(arrFiles[1], taReader)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root
Closing file, read 696635144 bytes in 163 transactions
Timings: Overall time = 106.982 s
Time per event: min=0.000 mean=0.009 max=6.802
As of right now, the dCache XRootD door is down - no - it’s back up
taDF <- galleryReader_df(taReader)
nrow(taDF)
[1] 20000
Set the scenario (we’ll need that again, so make it a function)
setScenario <- function(df) {
df %>% mutate(scenario = factor(fileEntry, levels=0:(length(scenarios)-1), labels=scenarios))
}
taDF %>% setScenario() %>% select(scenario, everything()) -> taDF
taDF
How many track action hits do we get per scenario.
taDF %>% group_by(scenario) %>% tally()
Let’s just get the death of the muons.
taDF %>% filter(status == 1) -> taDeathDF
We need the volume ID. I have a art FCL to run to get this information. I have files and they each may have a different set of volume IDs. Let’s try to do this in parallel. The examples in the help for parallel::mcparallel seem to be useful here.
Create the command strings,
runForVolIDString <- function(i) {
str_interp(
"PVS_CSVOUT=${csvout}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
list( csvout=scenarios[i], fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=arrFiles[i])
)
}
runForVolIDStrings <- sapply(1:length(arrFiles), runForVolIDString)
runForVolIDStrings
[1] "PVS_CSVOUT=everything_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root"
[2] "PVS_CSVOUT=noDipole_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089095_0/ARR_unified_noDipole_cyl.root"
[3] "PVS_CSVOUT=noInflector_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089097_0/ARR_unified_noInflector_cyl.root"
[4] "PVS_CSVOUT=noKickerQuads_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089099_0/ARR_unified_noKickerQuads_cyl.root"
[5] "PVS_CSVOUT=noKicker_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089123_0/ARR_unified_noKicker_cyl.root"
[6] "PVS_CSVOUT=noQuads_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089132_0/ARR_unified_noQuads_cyl.root"
Let’s run in parallel
jobs <- lapply(1:length(arrFiles), function(i) mcparallel(system( runForVolIDStrings[i], intern=T ), name=i))
res <- mccollect(jobs)
Now we need to load the csv files. Let’s do that in parallel too
jobs <- lapply(1:length(arrFiles), function(i) paste0(scenarios[i], "_volNames.csv") %>% read_csv(col_names=c("volumeUID", "volName")) %>% mcparallel(name=i) )
volNameTables <- mccollect(jobs)
names(volNameTables) <- scenarios
Let’s see what we got
volNameTables[["everything"]]
Join them together when the scenario name.
volNameTable <- map2_df(volNameTables, scenarios, function(df, sce) df %>% mutate(scenario = sce) )
volNameTable
Now let’s do lookups… [Note that I can’t do what I originally thought, which was split up the data frame into groups by scenario and then put in the factor by group – because in the end the rows get rejoined and then row-binded factor gets messed up.]
taDeathDF -> hold # in case we mess up
taDeathDF %>% inner_join(volNameTable) %>% select(scenario, eventEntry, volName, everything()) -> taDeathDF
Joining, by = c("scenario", "volName", "volumeUID")
taDeathDF
We can check this – we should see different volume names associated with volume IDs for different scenarios
taDeathDF %>% distinct(scenario, volumeUID, volName) %>% arrange(volumeUID)
Let’s plot where things die
taDeathDF %>% group_by(scenario, volName) %>% tally() %>% arrange(scenario, desc(n)) -> deathVolumeTable
deathVolumeTable
Did any of them actually decay?
taDeathDF %>% filter(volName == "ArcSection[00]")
Uh oh!!!! Hardly any decay!
taDeathDF %>% ggplot( aes(x = volName) ) + geom_bar() + theme(axis.text.x = element_text(angle=90, hjust=1))
Let’s look at the ‘everything’ scenario.
everythingFile <- which(scenarios == "everything")
I made a new TrackingAction reader
readr::read_file('trackingActionReaderCheckDecay.py') %>% cat
# Store the muon at death and whatever it decayed to at birth
from readGallery import GalleryReaderBase # Necessary for the base class
class TrackingActionWithDecayReader(GalleryReaderBase):
def __init__(self, inputTag, evLimit):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID',
'parentTrackID', 'volumeUID', 'status', 'x', 'y', 'z', 'px', 'py', 'pz', 'e', 'turn']
self.evLimit = evLimit
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))
def fill(self, ROOT, ev):
if ev.eventEntry() > self.evLimit:
return False
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2truth::TrackingActionArtRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2truth::TrackingActionArtRecord
for e in p: # Loop over elements and fill
# Let's only accept the muon at death or whatever the muon decays to at birth
if (e.trackID == 1) or (e.parentTrackID == 1 and e.status == 0):
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID, e.parentTrackID, e.volumeUID,
e.status, e.x, e.y, e.z, e.px, e.py, e.pz, e.e, e.turn ])
return True
Make an instance
taDecayClass <- createReaderClass_from_file('trackingActionReaderCheckDecay.py')$TrackingActionWithDecayReader
Load the data - let’s just look at one file since it takes a long time
taReader <- taDecayClass(artInputTag("artg4"), 20000)
getGalleryData(arrFiles[everythingFile], taReader)
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root
tadDF <- galleryReader_df(taReader)
nrow(tadDF)
[1] 52842
Make things easier to look at
tadDF %>%
mutate(r = sqrt(x*x+y*y),
pmag = sqrt(px*px + py*py + pz*pz) ) %>%
select(trackType, volumeUID, status, r, e, pmag, trackID, turn, x, y, z) -> tadDF
So I see a positron in alomst every case. But the muon appears to stop in the iron. Let’s just plot x,z
taDeathDF %>% ggplot(aes(x=z, y=x)) +
geom_point() + xlim(c(-10000, 10000))
tmp <- open3d()
ring <- cylinder3d( center=rbind(c(0,-90, 0), c(0,90,0)),
radius=7112,
sides=20, closed = F)
clear3d()
with(taDeathDF,
plot3d(x=x, y=y, z=z, type='p', main='Muon death positions',
xlab="x", ylab="y", zlab="z", xlim=c(-9000, 9000), zlim=c(-9000, 9000)))
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
This really doesn’t look like muons are being stored.
What fraction of muons have zero momentum when they die (so they die in something dense),
tadDF %>%
filter( trackID == 1, status == 1) %>%
summarize(pZeroPct = sum(pmag<0.1)/nrow(.)*100) -> pZeroPct
pZeroPct
So % of muons die in something dense. How many decay?
tadDF %>% filter(status == 1, trackID == 1, pmag>0.1) -> tad.muonsWithP.DF
tad.muonsWithP.DF
clear3d()
with(tad.muonsWithP.DF,
plot3d(x=x, y=y, z=z, type='p', main='Muon death with momentum positions',
xlab="x", ylab="y", zlab="z", xlim=c(-9000, 9000), zlim=c(-9000, 9000)))
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()
These really don’t look stored.
Look at the RingTrackingPlaneArtRecord.
readerClassSkel(extractClass = "gm2truth::RingTrackingPlaneArtRecord", fillClass = "gm2truth::BasicArtRecord",
writeFile = 'ringTrackingPlaneReader.py')
readr::read_file('ringTrackingPlaneReader.py') %>% cat
from readGallery import GalleryReaderBase # Necessary for the base class
class RingTrackingPlaneArtRecordReader(GalleryReaderBase):
def __init__(self, inputTag):
GalleryReaderBase.__init__(self, inputTag)
self.names = ['fileEntry', 'eventEntry', 'particle',
'trkid', 'parentid', 'pdgid', 'nturn', 'fracturn', 'detname',
'detnum', 'time', 'energy', 'x', 'y', 'z', 'px', 'py', 'pz',
'radius', 'angle']
# !!!! Modify the names accordingly - MUST MATCH self.vals.append CALL BELOW
def prepare(self, ROOT, ev):
GalleryReaderBase.prepare(self, ROOT, ev)
self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.RingTrackingPlaneArtRecord))
def fill(self, ROOT, ev):
validHandle = self.getValidHandle(self.inputTag) # Get the valid handle for gm2truth::RingTrackingPlaneArtRecord
if not validHandle.empty(): # Does it have data?
p = validHandle.product() # Get the corresponding data product (maybe a vector)
# Fill from gm2truth::BasicArtRecord
for e in p: # Loop over elements and fill
f = e.basicInfo
self.vals.append(
[ ev.fileEntry(), ev.eventEntry(), f.particle, f.trkid, f.parentid, f.pdgid, f.nturn, e.fturn, f.detname,
f.detnum, f.time, f.energy, f.pos.x(), f.pos.y(), f.pos.z(),
f.mom.x(), f.mom.y(), f.mom.z(), f.radius, f.angle ])
return True
rtpC <- createReaderClass_from_file('ringTrackingPlaneReader.py')$RingTrackingPlaneArtRecordReader
Make one and load the data
rtpR <- rtpC(artInputTag('artg4:RingTrackingPlanes'))
getGalleryData(arrFiles[everythingFile], rtpR) -> timings
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root
How many events did we see
length(timings$eventTimes)
[1] 10000
rtpDF <- galleryReader_df(rtpR)
nrow(rtpDF)
[1] 4
Only 4 events had a ring tracking plane hit?
rtpDF
Nathan pushed some now code that supposedly fixes this. Let’s look.
myNewFile <- '/home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root'
getGalleryData(myNewFile, rtpR) -> timings
Opening first file...
Successfully opened file /home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root
Closing file, read 765054 bytes in 57 transactions
Timings: Overall time = 0.130 s
Time per event: min=0.000 mean=0.000 max=0.016
length(timings$eventTimes)
[1] 1000
newRtpDF <- galleryReader_df(rtpR)
nrow(newRtpDF)
[1] 3203
newRtpDF
Now I see a lot more muons going around the ring!!
How many turns?
newRtpDF %>% ggplot(aes(x=nturn)) + geom_histogram(bins=100)
Not sure if this is right, but it looks better than before!
Let’s look at the volume names and see how mnay of these muons decay. It should be around 10% or so.
Let’s load the tracking action data and the volume name.
First, load the tracking action data…
tanReader <- trackingActionReaderClass(artInputTag("artg4"), 20000)
getGalleryData(myNewFile, tanReader)
Opening first file...
Successfully opened file /home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root
Closing file, read Timings: Overall time = 2.840 s
66328385 bytes in 597 transactions
Time per event: min=0.000 mean=0.003 max=0.685
tanDF <- galleryReader_df(tanReader)
tanDF
Just look at the muon death
tanDF %>% filter(trackID == 1, status == 1) -> tanDeathDF
tanDeathDF
Put in volume name
runVolString <- str_interp("PVS_CSVOUT=${csvout}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
list( csvout='nathanFixed', fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=myNewFile) )
runVolString
[1] "PVS_CSVOUT=nathanFixed_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 /home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root"
system(runVolString)
%MSG-i MF_INIT_OK: gm2 21-Mar-2017 14:31:10 UTC JobSetup
Messagelogger initialization complete.
%MSG
21-Mar-2017 14:31:11 UTC Initiating request to open input file "/home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root"
21-Mar-2017 14:31:12 UTC Opened input file "/home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root"
Begin processing the 1st record. run: 1 subRun: 0 event: 1 at 21-Mar-2017 14:31:12 UTC
21-Mar-2017 14:31:12 UTC Closed input file "/home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root"
TrigReport ---------- Event Summary ------------
TrigReport Events total = 1 passed = 1 failed = 0
TrigReport ------ Modules in End-Path: end_path ------------
TrigReport Trig Bit# Visited Passed Failed Error Name
TrigReport 0 0 1 1 0 0 pvsToFile
TimeReport ---------- Time Summary ---[sec]----
TimeReport CPU = 0.000000 Real = 0.000137
Art has completed and will exit with status 0.
volNames <- read_csv('nathanFixed_volNames.csv', col_names=c("volumeUID", "volName"))
Parsed with column specification:
cols(
volumeUID = col_integer(),
volName = col_character()
)
volNames
Let’s join!
tanDeathDF %>% inner_join(volNames) %>% select(eventEntry, volName, everything()) -> tanDeathDF
Joining, by = "volumeUID"
tanDeathDF
tanDeathDF %>% distinct(volumeUID, volName) %>% arrange(volumeUID)
tanDeathDF %>% group_by(volName) %>% tally() %>% arrange(desc(n))
So, about 1% decay!
tanDeathDF %>%
mutate(r = sqrt(x*x+z*z),
pmag = sqrt(px*px + py*py + pz*pz) ) %>%
select(eventEntry, trackType, volName, status, r, pmag, trackID, x, y, z) -> tanDeathdDF
tanDeathdDF %>% filter(volName == 'ArcSection[00]')
Do these correspond to the muons that make many turns? We need to do a join with the ring tracking planes
tanDeathdDF %>% inner_join(newRtpDF, by="eventEntry", suffix=c(".trk", "rtp")) -> tanDeathRTPDF
tanDeathRTPDF
Let’s look at guys that have more than one turn. Remember, we get a row for every hit. If we really want to see how many turns we get, we need to look at the last hit for the event.
tanDeathRTPDF %>% filter(fracturn >= 1) %>% arrange(eventEntry, fracturn) %>% group_by(eventEntry) %>% do(tail(., 1)) %>%
select(eventEntry, volName, r, pmag, fracturn )
Hmm - we’re missing some muons. Let’s look at those.
tanDeathRTPDF %>% filter(volName == 'ArcSection[00]') %>% arrange(eventEntry, fracturn) %>% group_by(eventEntry) %>% do(tail(., 1)) %>%
select(eventEntry, volName, r, pmag, fracturn )
What happened to event 791?
newRtpDF %>% filter(eventEntry == 791)
Let’s look at some muon tracks. How about 791 and 893. Since our file is on disk, this will be faster.